/* 
 * File:   fit.cpp
 * Author: mikdebian
 *
 * Created on July 10, 2013, 3:30 AM
 */
#include <stdio.h>
#include <limits>
#include <time.h>

#include <util/common/util.hpp>
#include <problems/boxcon/boxconbnb/lipbounder.hpp>
#include <problems/boxcon/boxconbnb/boxconbnb.hpp>
#include <problems/optlib/simpgrad.hpp>
#include <problems/optlib/gradboxdesc.hpp>
#include <problems/optlib/granradsearch.hpp>
#include <problems/optlib/locminchecker.hpp>

#include <lattice/lattice.hpp>
#include <tersoff/tersoffutils.hpp>
#include <lattice/diamondcubiclat.hpp>
#include <materials/silicon/silicontersoff.hpp>
#include <materials/silicon/siliconpowell.hpp>
#include <materials/carbon/carbonpowell.hpp>
#include <materials/carbon/carbontersoff.hpp>
#include "tersoffprob.hpp"

/**
 * Number of parameters to fit
 */
int N;
/**
 * Number of tries in Monte-Carlo
 */
int ns = 128;
/**
 * Treshold
 */
double trs = 0.001;

void makeBox(double* x, double eps, Box<double> box) {
    int n = box.mDim;
    for (int i = 0; i < n; i++) {
        box.mA[i] = x[i] - eps;
        box.mB[i] = x[i] + eps;
    }
}

class GradBoxDescStopper : public GradBoxDescent<double>::Stopper {
public:

    GradBoxDescStopper(double eps, int nstep) {
        mEps = eps;
        mMaxStep = nstep;
    }

    bool stopnow(double xdiff, double fdiff, double gnorm, double fval, int n) {
        bool rv;
        if (fval <= mEps) {
            rv = true;
        } else if (n > mMaxStep) {
            rv = true;
        } else
            rv = false;
        if (rv) {
            printf("found fv = %lf, nstep = %d\n", fval, n);
        }
        return rv;
    }

private:
    double mEps;
    int mMaxStep;
};

class GranRadSearchStopper : public GranRadSearch<double>::Stopper {
public:

    GranRadSearchStopper(double eps, int nstep, double mingrain) {
        mEps = eps;
        mMaxStep = nstep;
        mMinGrain = mingrain;
    }

    bool stopnow(int moves, double fval, double gran) {
        bool rv;
        if (fval <= mEps) {
            rv = true;
        } else if (moves > mMaxStep) {
            rv = true;
        } else if (gran < mMinGrain) {
            rv = true;
        } else
            rv = false;

        if (rv) {
            printf("found fv = %lf, nstep = %d, grain = %lf\n", fval, moves, gran);
        }

        return rv;
    }
private:
    double mEps;
    double mMinGrain;
    int mMaxStep;

};

/*
void fillInBox(Box<double>& box) {
    const double d1 = 0.9, d2 = 2.;
    TersoffParams tpar;
    SiliconPowell::fillInParams(tpar);
    double *p = (double*) &tpar;
    for (int i = 0; i < N; i++) {
        box.mA[i] = p[i] - BNBABS(p[i]) * d1;
        box.mB[i] = p[i] + BNBABS(p[i]) * d2;
    }
    printf("Initial values ==========\n");
    VecUtils::vecPrint(N, p);
    VecUtils::vecPrint(N, (double*) box.mA);
    VecUtils::vecPrint(N, (double*) box.mB);
    printf("=========================\n");
}
 */

void fillInBox(Box<double>& box) {
    TersoffParams tmin, tmax;

    tmin.mDe = 0.5;
    tmax.mDe = 10.0;

    tmin.mRe = 0.5;
    tmax.mRe = 5.0;

    tmin.mS = 0.5;
    tmax.mS = 5.0;

    tmin.mB = 0.5;
    tmax.mB = 5.0;

    tmin.mGamma = 0.5 * 0.0000001;
    tmax.mGamma = 3.0 * 0.000001;

    tmin.mNu = 0.1;
    tmax.mNu = 2.0;

    tmin.mC = 10000;
    tmax.mC = 200000;

    tmin.mD = 1;
    tmax.mD = 30;

    tmin.mH = -2;
    tmax.mH = -0.1;

    tmin.mLambda = 0.5;
    tmax.mLambda = 3;

    double *pmin = (double*) &tmin;
    double *pmax = (double*) &tmax;
    for (int i = 0; i < N; i++) {
        box.mA[i] = pmin[i];
        box.mB[i] = pmax[i];
    }
    printf("Initial values ==========\n");
    VecUtils::vecPrint(N, (double*) box.mA);
    VecUtils::vecPrint(N, (double*) box.mB);
    printf("=========================\n");

}

int main(int argc, char** argv) {
    if(argc != 4) {
        printf("usage: %s nparams ntries precision", argv[0]);
        exit(-1);
    } else {
        N = atoi(argv[1]);
        ns = atoi(argv[2]);
        trs = atof(argv[3]);
    }
    double eps = 0.00001;
    double refx[N];
    
    srand(time(NULL));

    TersoffParams tpar;

#if 1   
    // Silicon with Powell values
    // Fill in tersoff parameters with initial value proposed by Powell


    SiliconPowell::fillInParams(tpar);
    memcpy((void*) refx, (void*) &tpar, N * sizeof (double));

    // Fill in reference values
    TersoffRefVal refv;
    refv.mE = -4.630053;
    refv.mB = 0.975796;
    refv.mC = 0.511098;
    refv.mLC = 5.431;
    refv.mC44 = 0.797;
    refv.mKlein = 0.520;
#endif

#if 0  
    // Carbon with Powell values
    // Fill in tersoff parameters with initial value proposed by Powell

    CarbonPowell::fillInParams(tpar);
    memcpy((void*) refx, (void*) &tpar, N * sizeof (double));

    // Fill in reference values
    TersoffRefVal refv;
    refv.mE = -7.370;
    refv.mB = 4.408;
    refv.mC = 4.755;
    refv.mLC = 3.564;
    refv.mC44 = 5.778;
    refv.mKlein = 0.309;
#endif

#if 0  
    // Carbon with Tersoff values
    // Fill in tersoff parameters with initial value proposed by Powell

    CarbonTersoff::fillInParams(tpar);
    memcpy((void*) refx, (void*) &tpar, N * sizeof (double));

    // Fill in reference values
    TersoffRefVal refv;
    refv.mE = -7.422;
    refv.mB = 4.294;
    refv.mC = 4.819;
    refv.mLC = 3.561;
    refv.mC44 = 6.407;
    refv.mKlein = 0.229;
#endif

    TersoffParams tstore = tpar;

    // Constructing base lattice 
    DiamondCubic slbase(1, 3, refv.mLC);
    Lattice latbase;
    slbase.fillLattice(latbase);


    // Initialize utility classes
    TersoffUtils tutil(&tpar);
    tutil.setLatticeBase(&latbase);
    tutil.setEps(eps);

    // Check various components of the objective
    tutil.printValues();

    // Fill in bounding box
    Box<double> box(N);
    fillInBox(box);

    // Initialize tersoff objective (N first parameters to fit)
    TersoffObjective obj(&tutil, &refv, N);
    //    obj.getOptions() &= ~TersoffObjective::Options::INC_LATTICE_CONST;



    // Setup initial point
    //double xi[N], x0[N];
    //for (int i = 0; i < N; i++) {
    //    x0[i] = xi[i] = 0.6 * (box.mA[i] + box.mB[i]);
    //}
    // perform some scaling
    double aux[N];
    for (int i = 0; i < N; i++) {
        aux[i] = 0.5 * (box.mA[i] + box.mB[i]);
    }
    obj.setupScale(aux);
    obj.rescale((double*) box.mA);
    obj.rescale((double*) box.mB);
    printf("box a = ");
    VecUtils::vecPrint(N, (double*) box.mA);
    printf("box b = ");
    VecUtils::vecPrint(N, (double*) box.mB);
    obj.rescale(refx);

    //double v;
    //obj.rescale(xi);
    // VecUtils::vecPrint(10, (double*) xi);

    // Perform search
#if 0
    // Commented part
    TersoffProblem tprob(&obj, box);
    LipBounder<double> lb(&obj);
    lb.setLipConst(4 * N);
    BoxconBNB<double> bnb;
    bnb.addBounder(&lb);
    bnb.init(tprob);
    long long steps = 500;
    bnb.solve(steps);
    VecUtils::vecCopy(N, (double*) (bnb.getIncumbent().mX), xi);
    printf("v = %lf, steps = %ld\n", bnb.getIncumbent().mValue, steps);
#endif    

    double bestx[N];
    double bestv = std::numeric_limits<double>::max();
    for (int i = 0; i < ns; i++) {
        double xi[N], x0[N];

        for (int j = 0; j < N; j++) {
            double r = (double) rand() / (double) RAND_MAX;
            x0[j] = xi[j] = box.mA[j] + (box.mB[j] - box.mA[j]) * r;
        }
        double v;

#if 0
        // Initialize search procedure   
        GradBoxDescStopper stopper(0.0001, 400);
        GradBoxDescent<double> gbd(box, &stopper);
        gbd.setObjective(&obj);
        gbd.search(xi, &v);
#endif    

#if 1
        GranRadSearch<double> grs(box);
        grs.setObjective(&obj);
        GranRadSearchStopper stopper(trs * 0.5, 50000, 0.0000001);
        grs.setStopper(&stopper);
        grs.getProperties().mTresh = 0.1;
        grs.getProperties().mInitGrain = .1;
        grs.getProperties().mGrainDec = .1;

        grs.search(xi, &v);
#endif
       /*
        printf("%d: v = %lf, iadist = %lf, adist = %lf, iqdist = %lf, qdist = %lf\n",
                i, v, VecUtils::vecDistAbs(N, refx, x0), VecUtils::vecDistAbs(N, refx, xi),
                VecUtils::vecDist(N, refx, x0), VecUtils::vecDist(N, refx, xi));
         */
        Box<double> bx(N);
        double y[N];
        makeBox(xi, 0.1, bx);
        LocalMinChecker<double> lminchk(&obj);
        int nt = 1000;
        printf("Starting local minimum check\n");
        bool fl = lminchk.check(xi, y, nt, bx);
        if(fl) {
            printf("Local minimum confirmed\n");
            if(v <  trs) {
                obj.scale(xi);
                printf("v = %lf,  x = ", v);
                VecUtils::vecPrint(N, (double*) xi);
                obj.rescale(xi);
                printf("\n");
                fflush(stdout);
            }
        } else {
            printf("Fail to confirm local minimum.  Better value  %lf\n", obj.func(y));
        }

        if (v < bestv) {
            bestv = v;
            VecUtils::vecCopy(N, xi, bestx);
        }
    }
    //SiliconPowell::fillInParams(tpar);
    tpar = tstore;

    double kli;
    tutil.computeC44Base(kli);
    printf("Reference values: E = %lf, B = %lf, C = %lf, LC = %lf, C44 = %lf, kli = %lf\n vector: ", tutil.energyPerAtom(), tutil.computeBulk(), tutil.computeShear(), tutil.computeLatticeConst(),
            tutil.computeC44Base(kli), kli);
    VecUtils::vecPrint(N, (double*) &tpar);

#if 1
    /*
        memcpy((void*) &tpar, (void*) x0, N * sizeof (double));
        tutil.computeC44Base(kli);
        printf("Initial values: E = %lf, B = %lf, C = %lf, LC = %lf, C44 = %lf, kli = %lf\n vector: ", tutil.energyPerAtom(), tutil.computeBulk(),
                tutil.computeShear(), tutil.computeLatticeConst(), tutil.computeC44Base(kli), kli);
        obj.scale(x0);
        VecUtils::vecPrint(N, (double*) x0);
     */
    obj.scale(bestx);
    printf("v* = %lf, x = ", bestv);
    VecUtils::vecPrint(N, (double*) bestx);
    memcpy((void*) &tpar, (void*) bestx, N * sizeof (double));
    tutil.printValues();
    fflush(stdout);
#endif


    return 0;
}

